#例题2.8
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint,solve_ivp
from sympy import pi
def fun(t, w):
    x = w[0]
    y = w[1]
    dx = w[2]
    dy = w[3]
    # 求导以后[x,y,dx,dy]变为[dx,dy,d2x,d2y]
    # d2x为w[2]，d2y为w[5]
    return [dx,dy,-dy-3*x+np.cos(2*t),4*dx-3*y+np.sin(2*t)]
# 初始条件
y0 = [0,0,1/5,1/6]
yy = solve_ivp(fun, (0,100), y0, method='RK45',t_eval = np.arange(0,100,0.2) )
t = yy.t
data = yy.y
plt.figure(figsize=(12,6))
plt.subplot(1,2,1)
plt.plot(t, data[0, :])
plt.plot(t, data[1, :])
plt.legend(['x','y'])
plt.xlabel("时间s")
plt.subplot(1,2,2)
plt.plot(t, data[2, :])
plt.plot(t, data[3, :])
plt.legend(["x' ","y' "])
plt.xlabel("时间s")
plt.show()
